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Abstract 

. We propose a novel density based numerical method for uncertainty propagation under certain 

' partial differential equation dynamics. The main idea is to translate them into objects that we 

, call cellular probabilistic automata and to evolve the latter. The translation is achieved by state 

discretization as in set oriented numerics and the use of the locality concept from cellular automata 
theory. We develop the method at the example of initial value uncertainties under deterministic 
dynamics and prove a consistency result. As an application we discuss arsenate transportation and 
adsorption in drinking water pipes and compare our results to Monte Carlo computations. 



r^' ■ 1 Introduction 

> I 

^ The numerical treatment of differential equations that are subject to uncertain data has attracted a lot 

of interest lately. A prominent approach is to use Polynomial Chaos expansions [44j |8l [21] . It can be 
improved by decomposing the random space [42| . and only recently numerical implementations of this 
improvement have been investigated [T]. Alternative approaches are based on the Monte Carlo idea, 
J> i like the Markov Chain Monte Carlo Method [T^, Latin Hypercubc Sampling [32, the Quasi Monte 

'nI" ' Carlo Method [7] , importance sampling [5^ and the Multi-Level Monte Carlo Method [5] . Further well- 

known approaches use the Ito calculus [35l [23l [24l [37] or the Fokker-Planck equation [2^ . Although the 
approaches have proven to be successful for many tasks, they often encounter certain efficiency restrictions 
in higher dimensions of the random space. New methods are needed to meet these challenges. 
Time-continuous dynamical systems on continuous state space can be approximated by time-discrete 
. Markov chains on finite state space [T7]. This technique of state space discretization has led to the 

powerful tools of set oriented numerics [U [5] . It is especially useful to study ergodic theory, asymptotic 
dynamics and optimal control |28[ 111] . Recently, also contributions to uncertainty quantification have 
been made [T5] . 

In this paper we introduce a novel numerical scheme for uncertainty propagation in certain spatio- 
temporal processes. It is based on the concept of state space discretization and on ideas from cellular 
automata (CA) theory [501 [SI [13- We develop the method at the example of the propagation of initial 
value uncertainties under deterministic partial differential equation (PDE) dynamics and pave the way 
towards an extension to more general stochastic influences on the system. 

In particular we introduce a discretization of PDE which do not depend explicitly on the independent 
variables. First a finite difference scheme is applied to a PDE; the spatial and temporal continuum is 
replaced by discrete sites and discrete time steps. Second the state of the resulting system is discretized. 
Because this procedure emphasizes the interaction between neighboring sites, a property that strongly 
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resembles the locality and shift-invariance in cellular automata, the resulting completely discrete system 
is termed a cellular probabilistic automaton (CPA). Such an automaton is much simpler than the PDE 
and becomes accessible to very efficient simulation techniques. 

CPA basically consist of information about transition probabilities between discretized portions of phase 
space in a site's neighborhood. The transition probabilities arc interpreted as to approximate the evolution 
of the system's probability density in transfer operator theory [30] . Hence CPA can be used for uncertainty 
propagation [34] . While the translation from PDE into CPA may be rather time-consuming, the evolution 
of uncertainties with CPA is fast. The accuracy of the approximation depends on two parameters: one 
measures the state space resolution at every site, and the other the degree of locality, i.e. the extent to 
which correlations between neighboring sites are preserved. 

The paper is structured as follows. In Section [5] we formulate the problem of initial value uncertainty 
propagation under deterministic dynamics and deterministic boundary conditions. Here we also present 
the idea of density based uncertainty propagation through phase space discretization. By exploiting 
locality and shift-invariance of our problem this leads to the definition and discussion of CPA in Section 
[3j A consistency result for our construction is presented in Section lU In Section [5] we show how CPA 
can be extended to incorporate stochastic boundary conditions and apply the theory to the example 
of arsenate transportation and adsorption in water pipes. The results are compared to Monte Carlo 
computations. Finally we conclude in Section [6j 

2 Density Based Uncertainty Propagation 

In this section we first formulate the problem and then develop the idea of density based uncertainty 
propagation. Finally the CPA idea is derived in this context. 

2.1 Problem Formulation 

We are interested in the time evolution of uncertain initial data in a specific deterministic dynamical 
system. First we introduce some notation from probability theory and the Frobenius-Perron operator as 
the suitable tool to describe this process. Second we specify the deterministic dynamical system that we 
will work with, and third we formulate the problem. 

Let {X,A,ii) be a probability space, {X',A',ij,') a measure space and V : X X' a. random variable 
with distribution /iy. We say that V has density g if there is 5 e C^{X' , A' , fi') such that 



Let now {X',A',n') = (R™", A), where m,n e N, 6(R™") is the Borel cr-algebra and A the 

Lebesgue measure. A measurable map S : M™" — >■ K™" is called non-singular if X{S~'^{A')) = for all 
A' £ ;B(R™") with A(^') = 0. For any such map a unique operator can be defined on the basis of the 
Radon-Nikodym theorem [30] . 



Definition 2.1. Given a non-singular map S : M™" i-> R™", for g e £^(R'"") the Frobenius-Perron 
operator (FPO) Ps : £i(R"") ^ £i(R"") is defined by 



The FPO preserves positivity and normalization and hence describes how densities are mapped under 
phase space evolution. We focus on a particular type of phase space evolution. 




for all A' eA'. 



The set of densities on {X' , A' , ij,') is denoted by 



D{X') := {g e C\X',A\^i')\g > 0, |l.g|li = 1}. 
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Definition 2.2. Consider a deterministic dynamical system (T, R™", $) specified as follows: 

i) (T, +) is an additive half-group of time, 

ii) M""' = X/M" is the state space, where / = {1, m}, 

iii) the flow $ : T x R™" ^- R'"" is non-singular for all t G T, 

iv) there is a neighborhood U = {—r, ...,s} with r, s € No, r + s < m, such that $ has the locality 
property, i.e. that there is h : T x (R")l^l R" with 

V)i = h{t, Vi^r, Vi+s) 

for alH e T, w = (wi, Wm) £ X/R" and i G {1 + r, m - s}, 

v) and that the system acts as the identity on K ^ {1, r} U {m — s + 1, to}, i.e. = v\k 
for alH e T and all v e X/R". 

We will write := w) in the following and refer to [12] for further information on dynamical 

systems. Assume that there is a compact C R" such that il™ is positively invariant under the flow 
and fix T e T, T 7^ 0. 

Our main application is the analysis of a PDE 

df V = h{dxxV, dxV, v), v{x, t) e 

on a one-dimensional compact spatial domain x £ [a, 6] for a, 6 G R. Under certain regularity assump- 
tions a dynamical system like the above is obtained by applying a finite difference method with space 
discretization Ax = , where to G N,to > 2, and time step r. Then U is naturally induced by the 
choice of the finite difference scheme; e.g. usually U = {—1, 0, 1} is suitable to account for central second 
order difference quotients. Because of the PDE context we call / the set of sites. By only considering 
trajectories with v'^\k = k G x^fR", the system can be interpreted as to obey boundary conditions. 
The time evolution of uncertain initial data in the deterministic dynamical system is described by real 
random variables , ... : X ^ fJ™ on probability space {X, A, /i), where = We focus on 

deterministic boundary conditions: V°{x)\k = ke x i^R" for aU x e X. U V" has density e D{W^"), 
the density of V"'~^^ is given by application of the associated FPO: 5"^^ = ^$^((7"). The goal is to 
develop an algorithm that approximates the density evolution. It will be achieved by translating the 
system into a CPA in two steps. First the FPO is discretized via a state discretization procedure, and 
then locality and shift-invariance are used to further transform it into a CPA. 

2.2 State Space Discretization 

In this section first we introduce the concept of state space discretization. Second we investigate according 
densities, and third we construct a discretized version of the FPO. In principle these ideas are well-known 
in the literature [4j [5] . Here they are adapted to the special structure of the dynamical system. 

Definition 2.3. A partition or coding of is a finite collection of disjoint sets {fieleeE whose union 
is n. We call e E E the symbol of coding domain He, and the coding map is the function T : H> i?, 
where T{v) = e if u G ^le- A partition is called uniform, if there is a resolution Afi S R such that fig is 
a n-dimensional hypercube with side length Afi for all e € E. 

To avoid technical complications in the following proofs we only consider uniform partitions while devel- 
oping the theory. They are also the ones that are relevant in practical algorithms. 

A partition _E of 17 with coding map T and \E\ ^ N naturally induces a partition E^ of f]™ with coding 
map 

f : 17" ^ E',v^ f{v) with {f{v)\ = T{vi) for i e /. 
Note that \E^\ ~ mN . For Lp G E'\ where J C /, we write 

fi^ = {z;efi™|Vje J: f{v)U)^v{j)}. 
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Now we study densities that are compatible with state space discretization. For this purpose we introduce 
the measure space {E^ ,'P{E^),j), where V{E^) is the power set of E^ and 7 is the counting measure. 
The densities D{E^) consist of the weight functions 

5 : ^ [0,00], 9{f)=P^, 

where (pip)ipeB^ ^^rc nonnegative numbers with '^^p^E' Pv ^ 
Definition 2.4. 

i) £^(M'"") = span(i3) is the finite-dimensional £^ (M"'")-subspace of piccewise constant functions 
with basis B = {xn,^ / ^{^if)},pi£E' ■ The set of piccewise constant densities is given by Z3^(E'"") 

ii) The coordinate representation kb '■ £^(M"'") g i-> with respect to the basis B is 
given by ^sigK'p) = for ^ e E' and g = Y.i,eE' W^o^- Obviously kb(D^(K'"")) = D{E'). 

iii) Let p € E^ such that pi = T{ki) for all i £ K. The densities that are compatible with the boundary 
conditions are given by 

Dbc{E') := {g € D{E') \ g{^) = if ^ p}- 

By averaging in the coding domains every function in L^iW^"") can be mapped to a piecewise constant 
function. 

Definition 2.5. A restriction operator to the subspace of piecewise constant functions is given by 

R : /:i(M'"") ^ 4(M'""), R{g) = ^ TT^Xa,, 

where 

g{w)dw. 



R is idempotent, i.e. R o R ^ R, and furthermore i?(£)j,(R'"")) C Df{R""^). In the following we wiU 
use the restriction operator to construct a discretized version of the FPO on density level: RPg,T . This 
procedure is well-known in ergodicity theory when invariant measures are approximated. There it is 
called Ulam's method [40] . 

The matrix representation of the linear i?P$T j^i (]g™„-) is given by Pb = kbRP<!>^ e with 
entries 



PB,ip,ii) is the probability of finding a realization of a random variable with uniform density in VL^p in 17^, 
when <I>'^ is applied. Hence we may interpret Pb,lp,^ as the transition rate from D.^p to of a finite 
state Markov chain on {^if,}^p£E' ■ This chain approximates the behavior of the dynamical system for 
uncertain initial values. 

In the following we regard Pb ■ Dbc{E^) — >■ Dbc{E^) as a function which maps densities that are 
compatible with the boundary conditions by matrix multiplication. 



2.3 Using Locality - Towards Cellular Probabilistic Automata 

E^ grows exponentially in m. For a growing number of sites it becomes numerically expensive to obtain 
global transition rates and to handle global states and densities. 

However, our dynamical system has a special structure: We use the locality property to approximate 
the set of global transition probabilities by several identical sets of local ones. This is possible out of 
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two reasons. First because we find identical dynamics at all sites away from the boundaries, and second 
because the transition probabilities at one particular site mainly depend on the state of its neighborhood 
rather than on the whole global configuration. 



For the formal definition of these local transition probabilities we need to introduce the shift by / G Z on 
finite grid J C Z. It is given by 

where F is an arbitrary set, e. g. F = E or F = D{E^). Moreover, for arbitrary V = {—p, ■■■,q}, W = 
{—t, u} with p, q,t,u E No and I E Z wc use the conventions I + V = {—p + I, ...,q + 1} and V + W ^ 
{-p-t,...,q + u}. 

Definition 2.6. Let V = {— p, q} with ^,^£1^0 and p + q + r + s<m. A local function /o : E^~^^ — > 
D{E^) is then given by 

hmw) = TTT^ ^ , 

where i = l+ p + r,ipe E^^^ and e 

Note that because of the locality property the definition is independent of the chosen site i G { 1 + p + 
r, TO — q — s}. The set V controls the degree of locality, i.e. the number of sites that give rise to a local 
transition. It will turn out that by enlarging it we can diminish the error of the locality approximation. 
In the following section we develop a method of how to combine several such local transitions to approx- 
imate a global one. This will finish the construction of a CPA from the FPO. 



3 Cellular Probabilistic Automata 

CPA are defined by extending the definition of deterministic CA according to [6l [16] : In CPA the local 
transition function specifies a time- and space-independent probability distribution of next states for each 
possible neighborhood configuration. Because we do not want to follow one realization but rather the 
whole ensemble, unlike in the literature we define CPA to work on densities. This enables their utilization 
for uncertainty propagation. 

In the last chapter we showed how the discretized FPO Pb on state space Dbc{E^) can be used to 
approximate the FPO P$t on _D(M'""). CPA further approximate the discretized FPO on a product 
space of local densities, see Fig. [T] for a sketch. Uncertainty propagation with CPA therefore requires 
two definitions. First one about how to translate between global densities and the product space of local 
densities, and second one about how to evolve local densities in time with the help of the local function. 
Because the definitions can be best understood for V = {0}, in Section [01 we first introduce CPA in 
this special case to demonstrate the basic construction. Afterwards we develop the de Bruijn calculus as 
a connection between local and global objects for more general V in Section [321 This connection leads 
to the generalization of CPA to general V in Section 13.31 



3.1 Cellular Probabilistic Automata: A Special Case 

A crucial step is to translate between global densities Dbc{E^) and a (subset of a) collection of local 
densities {D{E^))^ , where / contains the sites away from the boundary. We introduce an operator 
/§ : Dbc{E^) {D{E^)y and a concatenation operator a : {D{E^)Y Dbc{E^)- P localizes 
the information to densities on states of length V and thus erases far-reaching correlations, a in turn 
constructs global densities out of information about local densities. As we will see in the next section, 
this process is by no means unique and requires some technical refinement of the space of local densities. 
However, for V ~ {0} there are canonical definitions for a and (3: multiplication of local probabilities for 
independent events and calculation of marginal distributions. 
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Figure 1: The relations between the FPO P$t and its approximations. By state discretization we obtain 
the discretized FPO Pb which stiU works globahy, and by exploiting locality we approximate Pb further 
by the CPA with global function /. The state space on which the CPA operates is a collection of local 
densities, see text. 

Definition 3.1. As before let / = {1, to}, U = {—r, s} for r,s,m £ Nq with 1<to, f + r + s<TO, 
and K = {l,...,r}U {to — s + 1, m} . Furthermore / = {1 + r, ...,m — s}. 

i) We set a : {D{E)y Dbc{E'), g a{g) with 

ii) and /3 : Dbc{E') -> {D{E)y , g ^ j3{g) with 

XeE' s. t. x{i)=e 

Definition 3.2. A cellular probabilistic automaton (CPA) is a tuple (/, U, E, /o), where for m,r,s £ No 
with 1 < m and f + r + s < m 

i) / = {1,...,to} is a finite grid, 

ii) U — {— r, s} is the neighborhood, 

iii) i? is a finite set of local states 

iv) and /o : E^ — > D{E) is the local function. 

With the boundary conditions p £ E^ the global function is given by 

f : {D{E)y ^ {D{E)y , g^fig), 

figmw = E n 9ij)Mj))foi^m- 

^eE" s.t. y,{h-i-) = p{k) j&(i+U)\K 

for k£K(M+U 

The trajectory starting with G {D{E))^ is given by the sequence (g")„gN, where = /(.g"~^) for 
n G N+. 

A CPA can be used to evolve an input distribution f3{g) for g G Dbc{E^) via the global function. 
After n time steps the approximated global density is then given by af"'(3{g), see also Fig. [1] with 
= D^^g = {D{E)y. The role model for the global funct ion is the matrix operation with the 
discretized FPO Pb'- The product of the transition probability with the probability of being in a preimage 
state is summed up over all possible preimage states. A probability is assigned to a preimage state G E^ 
by multiplication of local probabilities like in the definition of a. 

To cope with boundary conditions it is necessary that the global function only operates on / instead of 
/. Deterministic CA are special cases of CPA: assume that for all (f G E^ there is e G such that 
.fo{^){e) = 1 and that the input is deterministic. 



3.2 De Bruijn Calculus 

To generalize the construction to arbitrary V wc first study the relation between local and global objects 
in more depth. We introduce the de Bruijn density calculus on the basis of pattern ideas in cellular 
automata theory |14 1 141 j . in the theory of de Bruijn graphs |39j and in pair approximation |22) . 
As before, we introduce an operator that localizes the global information to densities on states of 
length V, this time \V\ > 1. The precise definition of is still rather straight forward: marginal 
distributions dismiss all information but the one over a certain range V. We will find below that the 
reconstruction of global densities out of local information by aw is more involving. However, let us first 
define /3w- 

Definition 3.3. Let V = {—p,...,q} and W — {—t,...,u} for p,q ^ N and t,u E Z, — t < u. : 

D{E^+^) -> {D{E^))^ is given by 



Example 3.1. het E = V ^ W = 1}, p G (0, 1) and g,~g e D{E^+^) be given by 

( p if V (001) 

5(^) = < 1-P if V- (100) , 5(^) = 

[ else 

We find that I3w{g) ~ l^wig) with 



p{l 


-p) 


if ?A = 


(000), V' = (101) 






if = 


(001) 


(1- 


-p? 


if = 


(100) 







else 





P 


if = 


(00) 


r (i-rt 


if V = 


(00) 


1-p 


if = 


(10) , 


Pw{g)iim = { P 


if ?A = 


(01) 





else 




[ 


else 





Ex. 13.11 shows that information is lost under I3w, i-e. different global densities are mapped to the same 
collection of local densities. Now we are interested in aw '■ D{E^^^) — )- {D{E^))^ . Although the 
properties of /3w allow to define aw as the solution of a linear nonnegative least squares problem |31| . 
this algebraic approach is not appropriate. We rather suggest a probabilistic approach that fulfills two 
requirements: first aw shall degenerate to simple multiplication of local densities for V = {0}, and second 
aw and (3w shall be inverse on important sets. For this purpose we first introduce several definitions, 
see also Fig. [2al 

Definition 3.4. 

i) = (r{E^))^ is the set of de Bruijn states, where r{E^) is the power set of E^ . The elements 

are called patterns of size V. 

ii) The subset of extendable de Bruijn states is given by 

xYbc = {* e XdB\y^ e wyip e $(i)3Vj e E^+'^yi e W a,(V')k = (^}, 

the set of completely extendable de Bruijn states, 
ni) D% = {D{E^))^ is called the set of de Bruijn densities, 
iv) The subset of extendable de Bruijn densities is given by 

DYbc ^{g^Djsl ^^ew supp g{i) £ XYbJ- 

The idea behind extendable de Bruijn states is that every pattern can be extended to a global state by 
gluing suitable patterns on it. An example of an extendable de Bruijn density is /Swig) in Example 13.11 
for example pattern (10) at site can be extended by (01) at site 1 to the global state (101), because the 
patterns coincide in the overlapping state 0. We find that this observation can be generalized. 
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Figure 2: (a) The relation between global densities and de Bruijn densities, (b) An example of an 
approximation for W = {—1, 0, 1}, V — {0, 1} and i = 0. The thin and dark grey box that covers sites 
and 1 is the factor at site 0. The medium box that covers sites —1 to 2 is the factor at site 1: the 
local state at site 2 depends on the local states from —1 to 1 (medium grey part). In the approximation 
(dashed line) the local state at 2 only depends on the one at 1. The thick box covering sites —1 to 1 is 
the factor at site —1. The local state at site —1 depends on the local states on sites and 1 (light grey 
part). In the approximation the dependence stops again at the dashed line. 



Lemma 3.1. im (fiw) ^ D^Be- 

Proof Let g € im {/3w), j G W and (p e with g{j){'p) > 0. Then there is g G D{E^+^) and 
tp £ E^+w such that g(V') > and crj(V')|v = ^- But furthermore already g{i){cri{il})\v) > for all 
i e W, and therefore Xigv^supp g{i) G X^^. □ 
Because only the extendable de Bruijn densities are addressed by Pw, we only look for aw on them. 
Our choice is motivated by the following calculation, for which we first introduce some notation. For 
J C Z and A,B (- E"^ , fj,[A] = J2<peA9i'i^) denotes the distribution associated with g £ D{E'^), and 
/.t[A|_B] = ^^^g^^ the conditional probability. Furthermore 

{^\j} = {^eEJ\p\j = ^\j} 

for J, J C Z, J C J, J C J, and t/j G E\ and we also write {"01 j} = l^"! j| if J is clear from the context. 
We calculate iovieW and ip G E^+^ that 

mm =/i[{?M{-t i}+v}] 11 



M[{V'l{-t,...,z}+y_}] 

u 

n l^l{^\{k-p}}\M{K..a} + V+}]^^[{iJ\^+v}] Yl MKV'I } I {V'l {-t,... J}+V- }] . 



i-1 



k=-t l=i+l 

where = {—p + 1, q} and V- = {—p, q — 1}. If there are no far-reaching correlations, we expect 
the following approximations to be suitable: 

MlWkfc-p}} I M{k,...,^}+v+}] ~ MKV'kfc-p}} I {V^lfc+y+j], 

M[{V'l{i+<?}} I {i^\{-t,...j}+v-}] ~ ^^[{^\{l+q}} I 

for k G {— t, j — 1} and I G {i + 1, u}. They resemble a Markov property of order — 1 in space, 
see also Fig. [2bl a site's state is independent from the states at sites that are more than \V\ — 1 sites 
apart. 
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Lemma 3.2. Let fij be the distribution associated with I3w{9){j) S D(E^) for j e W. Then 
fi[{ip\{k-p}} I = M/c[{CTfc(V^)|{-p}} I {'7fe(V')|v+}], 

^^[{^\{l-p}} I {V-li+y-}] = ^^l[WlW\{^p}} I {o-iWk-}], 

for i e VK, fc e {-t, i - 1} and I e + li}. 

froo/ Without loss of generahty we prove the statement only for k e {— i — 1}: 

A^ii^lffc-p}! I {V\k+V+i\ = = ^ 



Using the approximations introduced above and the lemma we find that 

i — l u 

■i— 1 u 
k=-t l=i+l 

This way we are led to the following definition. 



Definition 3.5. Let i e W. Then aw^t ■ D(E^+^) is given by 



i-i 



where /Xj is the distribution associated with g{i) S D{E^) for j e W , and -0 S 

For W = {u}, u E Ij, the definition simplifies to (5) (V') = 3(^t)(o'u(V'))- For = {0} the conditions 

vanish, and we get back simple multiplication, awAig){'>P) — IlfceH' We justify the general 
definition in the following lemma. 

Lemma 3.3. Let g G and i £ W. Then awAg) & D{E^+^). 

Proof Let g E -D^^ and i E W. aw,i{g) ^ because there is at least one extendable pattern with 
nonzero probability. We prove that it is also normalized in the following. Without loss of generality we 
assume that i = u. Then 

u-l 

XI OtW,u{9){'f) = XI W ^J■k{{(^k{'f)\{~v}}\{^k{'p)\v+}\^J'u[{(Tu{'p)\v}] 
^(zEV+w ^(zEV+w k=-t 

= X 5Z M-t[{CT-i(<p)|{-p}} I {cr-t(<^)|y+}] 

u-l 

n ^^k[Wk{(p)\{-p}}\{(Jk{'f>)\v+}]^J'u[{(yu{(p)\v}\ 
k=-t+l 
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u-1 

= n ^^k[{crki'f)\{-p}}\{<Jk{'f)\v+}]^^u[{<Jui^)\v}] 

0g£;V + H.'\{-t} k=-t + l 

The steps indicated by ... follow by induction in □ 
We find the following results for our construction. 

Lemma 3.4. Let i G W and g G im [fiw]- Then aw,i{g) = ctwjid) for all i,j G W. 

Proof We show that aw,i{9) = ctw,i+i{g) for all i £ W\{—t}. An index shift to the left can be proven 
analogously. 

Note that v}] ~ and that for k e {— <, « — 1} and I E {i + 1, u} 

rr /,M II r MM 11 .9(fc) K (V') I V ) 



MWim{,}}\Wimv.}]^ 



Therefore avy_i(g)('0) and aw,i+i{g)i'<P) have the same numerator and only differ in the denominator. It 
is enough to show that a factor in the denominator may be shifted one step to the right: Let i E W\{u} 
and g = Pwig) for g € D{E^+^). Then 

E ^(y') = E 

= E E = E /3>^(5)(* + l)(x) □ 

Theorem 3.5. Let g e im(/3vi'). Then I3wcitw,i{,g) = .9 for all i E W. 

Proof Let i E W and g € im{/3w)- We prove = g{j) without loss of generality only for 

j = u. 

With Lm. 13.41 and because a conditional distribution is a distribution as well, for -0 E 



Pwaw.i{9){u){ij) =Pwaw.u{9){u){ij) = E Oiw,u{9){V') 

= E E M-t[{f^-t('^)i{-p}} I {f^-t('^)k+}] 



n /^A;[{fTfe((^)|{_p}} |{crA;(<^)|y_,}]Af„[{cr„((^)|v}] 

E n ^^k[{(Tk{'f>)\{-p}}\{(Tki'f)\v+}]^^u[{(Tu{^)\v}] 

9{u){ip)- 



fe=-t+i 
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The last steps follow by induction in □ 
Now we focus on global densities that are preserved under aw,il3w for all i £ W. We will see that they 
enable an algebraic interpretation of 13^. Our choice of a^A is further justified thereby. 

Definition 3.6. g e D{E^+^) is called T^-factorizable, if g = aw.iPwig) for aU i G W. 

An example of a {0, l}-factorizable density is g in Ex. 13.11 g in the same example however has correlations 
over 2 sites and is not {0, l}-factorizable: a state only has positive probability if the local states at sites 
and 2 differ. The correlations are ignored by /3w- 

Lemma 3.6. Let i E W and g G im (fiw)- Then aw,i{g) is V-factorizable. 
Proof Let g e im {/3w) and i,j eW. Then by Lm. [SHand Thm. [33] 

aw,i(5) = o^wjig) = oiw,jPwoiw,i[g)^ 

and hence aw,i{g) is V-factorizable. □ 

Tiieorem 3.7. For aU g e D{E^+^) there is a TZ-factorizable g <E D{E^+^) such that pwig) = Pw{g)- 

Proof Let g G D{E^^^) and choose g = aw.uPw{g)- Then g is y-factorizable by Lm. 13. 6[ and the 
definition does not depend on the site u. Furthermore /3vy (5) = Pw(^w,uPw{g) = Pw{g) by Thm. 13.51 □ 
Pw induces equivalence classes on D{E^^^) by collecting all global densities with the same image in one 
class. According to Thm. 13. 71 each equivalence class contains at least one y-factorizable state. Moreover, 
we know that there is exactly one such state, because j3w is injective on these states by definition. It 
is given as the image under aw,iPw of any state in the class for any i G W . Therefore it is possible to 
choose the F-factorizable states as the representatives of the equivalence classes. These representatives 
are preserved under aw,iPw for any * S W . Ex. 13.11 provides an example: g and g are in the same 
equivalence class, and g = awfil^wig) is the unique F-factorizable representative of the class. 
However, in general aw,i[g) is not F-factorizable if g G £'2Be\im {Pw)- There is a degree of freedom in 
how to map a density collection to a global density on this set. We choose the arithmetic mean over all 
aw,i, where i G W . Note that for g G im {Pw) the definition then coincides with any aw,i- 

Definition 3.7. aw ■ -> D{E'^+'^) is given by aw{g) = T.tt^w ^w.i{9)- 

It is clear that aw{g) is a density by similar reasoning as for awsig)- 

3.3 General Cellular Probabilistic Automata 

With the de Bruijn calculus at hand we can now generalize the definition of CPA to general V . 

Definition 3.8. As before let / = {1, m}, U ~ {— s}, V = {—p, q} for p, q, r,s,m E Nq with 
m > 1,1 + p + q + r + s < m and K = { 1 , . . . , r} U {m — s + 1 , . . . , m} . We now set ii = I + p + r, 
ir = m — q — s and / = {i/, v}. 

i) We set a : 1?^^^ Dbc{E^), 9 ^ a{g) with 

a(Q)(ib) = i "/(£')(V'l{l+r,...,m-s}) if i^\K p 

I else 

ii) and /3 : Dbc{E') ^ i^Le, .9 ^ Pig) with 

kgmW = E aix)- 

x\i+v=o--i{i') 

Definition 3.9. A cellular probabilistic automaton (CPA) is a tuple {I,U,V, E, fo), where for 
TO, p, q,r, s E No with to > 1 and l+p + q + r + s<m 
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/I 23456 78 
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Figure 3: An example of a CPA with / = {!,..., 8}, U = {-1,0,1} and V = {-2,..., 2}. The global 
function considers patterns located at / = {4, 5} and sketched in the image (lower part) by rectangles of 
small and large height, respectively. The corresponding preimage patterns (upper part) are larger due to 
the neighborhood, and their probability of occurence is influenced by boundary conditions (dashed parts). 
From an implementational point of view, V may be constructed from V = {— 2,...,1} and W = {0, 1}, 
see Section 15.21 Focus on a pattern at site 4. The transition probability from a preimage pattern is 
calculated from the information about two subtransitions between light and dark grey subpatterns. 

i) / = {1, m} is a finite grid, 

ii) U — {—r, ...,s} is the neighborhood, 

iii) V = {— p, g} gives rise to de Bruijn patterns, 

iv) E is a finite set of local states 

v) and /o : E^+'^ D{E^) is the local function. 

With the boundary conditions p € for K = {1, r}U{m— s+1, m} and the definitions ii = 1+p+r, 
ir ~ m — q — I ~ {ii, v}, U (i) ~ {—"(*): and u,u : I ^ U with 

( i - ii if I e {ii, ii+r -1} 
\ r else ' 

( ir — i Hie {ir — S + 1, ir} 
\ s else ' 

the global function is given by 

f ■■ D^Be ^ D'jBe, 9 ^ I [q) , 

f{9)ii){i')= Oiu{i)i'^d9)\u{t))i'P\v+u{z)) ■ foi'P)W- 

for keKni+U+V 

The trajectory starting with G {D{E^))' is given by the sequence ((7")„gN, where g" = f{g"~^) for 
n e N+. 

See Fig. |3]for a sketch of how the CPA works on general patterns. Remark that /o is not arbitrary but 
connected to a dynamical system with the locality property. By investigating this relation we can assure 
that the global function is well-defined. 

Lemma 3.8. fiD^sJ C D',^^. 

Proof Because D^^^^ ~ {D{E^)y for \V\ — 1, the statement is trivial in this case. So we focus on 
\V\ > 1. We show that without loss of generality any pattern in the support of any site in the image can 
be extended to the right by a pattern in the support of the neighboring site. 

Let g S -D^Be'* S {"ill ■■■■,'ir — 1} and ip € supp f{g){i)- By the construction of / and a^/jj-j we know 
that there is <y9 G E^^^ such that <7j{{p)\v £ supp(5)(i + j) for all j G U and fo{tp){ip) > 0. Because 
of the extension property of the preimage g we can find (p G E^^^ such that aj{(f)\v = aj+i{cp)\v G 
supp {g){i + 1 + j) for aU j G {— r, s — 1} and as{^p)\v & supp {g){i + 1 + s). 



u{i) = 



12 



This enables us to find G supp f{g){i + 1) such that V'lv'- = o'i('0|v+)) as we will show in the following. 
Hence the pattern ip may be extended to the right by "0, and the proof is complete. 
Since fo{^p){il') > a-nd the partition is uniform, there is an e-ball with respect to the 2-norm in R"'", 
e > 0, such that 

Because the set is just restricted on sites i + U + V due to the locality property, we may independently 
restrict at site i + 1 + q + s and can still find e' > with 

In the second line we have again used the locality property, and the equality sign holds due to 
(^-iifW^+u) = o'-(i+i)(<^|v_+£/)- Wc now define V G by ip\v^ — (Ji{4>\v+) and choose '4i{q) such 
that there is e" > with 

Therefore 

/o((^)(Vi)>0, and V^G supp + ^ ^ □ 

We denote the case of ii = v with Fmax and find that al3{g) = g for all g E Dbc{E^) and $a{g) = g 
for all g € {D{E^)Y . Furthermore ~ {0} in this case, and it can be calculated for g E D^dBe and 
ip e E'^+^ that 

a{o}{<J^,{g)){^p\v) = g{ii){^\v)- 

For Ip € E"^ then 

where the sum is taken over all ip E E^^^ such that the boundary conditions are fulfilled. For Vmax the 
evolution of the global density is calculated directly, and locality is completely omitted. 



4 Consistency 

The goal is to prove in Section r4.3l that time evolution of probability densities under CPA can approximate 
time evolution under FPO arbitrarily close. We prepare this consistency result by investigating first the 
two features state space discretization and locality in more depth in Sections 14.11 and 14.21 respectively. 



4.1 State Space Discretization 

In this section we just compare the discrctizcd FPO to the real FPO without considering locality. There 
are many ways to study distances between probability measures [5] . In our density based formulation we 
use the £^-norm for probability densities which leads to the notion of strong convergence in the literature 
|30j . It is well-known that in this norm the discretized operator converges pointwise to the FPO for 
increasing state space resolution in the case of Lipschitz continuous input densities |26| . But because 
the image under RP^t is in general not continuous, for iteration we need to generalize the result to 
£"'^-functions. We start with some necessary tools. 

A subset M C X of a metric space {X, d) is called totally bounded, if for every e > there exist n G N 
and xi, ...,Xn £ M such that 

n 

M C |J{a; e X\d{x,x{) < e}. 

i=l 
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Totally bounded subsets of £''(M") can be alternatively characterized in a functional analytical sense by 
the theorem of Kolmogorov-Riesz [151 145]. In the following |.| denotes the 2-norm in M". 

Theorem 4.1. Kolmogorov-Riesz 

A set M C £^(]R"), 1 < p < oo is totally bounded if and only if the following criteria are fulfilled: 

i) M is bounded in £^'(M"), i.e. sup^^j^j \\g\\p < oo. 

ii) For every e > there is some R so that for every g & M 



\giv)\Pdv < eP 



\v\>B. 



iii) For every e > there is 5 > so that for every g E M and w E M" with |w| < 5 

\g{v + w)- g{v)\Pdv < e^. 



Theorem 4.2. Let g e Z)(M™") with supp(.g) C fi™ and T : ft E a uniform partition with resolution 
Ari. Then R converges pointwise to the identity with respect to the £^-norm, 



l-R(.9)-.9lli ^0 (AO^O). 



Proof 



\Ri9)-g\\ 



f Xn,(«)^7^ / g{w)dw- Xn^iv 



)9iv) 



<f>eE 
1 

AO™" 



dv 



<y — 
-y — 



\g{w) — g{v)\ dwdv 



\g{v + u) — g{v)\ dudv 



The last step involves a change of variables from w to u :~ w ~ v, and flip — v :~ {u — v \ u E ilip}. As 
for V e il^ with ip G , we calculate with Fubini's theorem 



R{9)-9\\i < y / / 



< 



V£E 
1 



^lG[-A^,A^^]'' 



\g{v -\- u) — g{v) \ dudv 



Af7" 



tie [-An, An]™" Jug 



\g{v + u) ~ g{v)\ dvdu. 



Let e > 0. Because the set {g} C £"'^(]R") is totally bounded, the theorem of Kolmogorov-Riesz, Thm. 
4.1[ guarantees that there is (5 > such that 



/ \9{v + u)~g{v)\dv<-^ 

J-uGR"" ^ 



/■uG 

for ImI < 6. If we choose Afl such that AD, < / we ensure that 



\ yM'< 



Y = VrnnAn < S 
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for all ue[-An, AfJ]™" and hence 

\\R(.9) - 9\\i < ir^ — / -^du = e. 

Therefore - g|l i ^ for Af} ^ 0. □ 

4.2 Locality 

In this section we investigate in more depth the role of locality in approximating the discretized FPO Pb 
by a CPA. It turns out that the CPA covers the dynamics of the underlying if only the support is 
considered. However, we cannot obtain the precise behavior of the global density with CPA in general, 
as will be shown with two examples. In the end we will see that according errors vanish for maximal 
pattern size. 

Let us start with the cover property. 

Lemma 4.3. For all n e N, for all g e Dbc{E^) and all i e I it holds that 

supp WPSigm) C supp (r/3(g)(*)). 

Proof Let e N, g E Dbc{E^), i G / and x e supp {j3P'^{g){i)) e E"^ . Then there is V e -E^ 
such that PB{g){ip) > and ai{'ip)\v — X- Let ipn ~ ip. Per induction it can be shown that we 

can find ipQ, (pn—i G such that PB,ipk-i,ifik ~ ^T^TT ) — 

^ > and > for 

k E {1, ...,n}. Because for all j E I 

we conclude that fo{(Tj{(pk-i)\u+v){<^j{'Pk)\v) > for aU j E I. Furthermore f3{g){j){aj{ipo j) > for all 
j E /, and therefore f P{g)ij)icrj{(pi)\v) > for all j E I. This induces P K9){j){'^j{'^2)\v) > for all 
j E I and so on, and therefore ${g){j){crj{(pn)\v) > for aU j E I. RecaUing that cri{fn)\v = Xi 
conclude that x G supp {f^${g){i)). □ 
However, we cannot recover the precise global behavior of the discretized FPO from a CPA. The errors 
that can occur are twofold, and we will provide examples for both types here. On the one hand it may 
happen that correlations over \V\ sites are not preserved because we work on patterns of size V. On the 
other hand we will see that even ioi U C V in general there are locally allowed transitions of a global 
state that are not allowed in a global consideration with Pb- This is remarkable, since such behavior 
was ruled out for the underlying dynamical system by the locality property. But because a CPA only 
computes locally, that may also lead to errors. While the first error type is a true locality effect, the 
second arises from the interplay of locality and state space discretization. 

Example 4.1. This example shows that in general correlations over \V\ sites are not preserved. We 
compare one CPA time step to one time step with the discretized FPO. Consider / = {1, 2, 3}, and the 
dynamical system that is given by the identity on ft™ ~ [0, 1]^, i.e. U ~ {0} and K = 0. We choose 
the partition Vlo = [0,0.5) and = [0.5, 1], i.e. E = {0, 1}, and look at the CPA with V = {0, 1} and 
/ = {1,2}. We find /o(<^)(V') = ^v-i' ^^r all e E^"'^ . 

Consider g E Dbc{E') = D{E^+^) with W ^1 from Ex. |3lll /3(.g) = Pjig), and also //3(g) = I3j{g). 
So a${g) = aj(3j{g) = g. However, Psig) ~ g, and so af${g) ^ Psig)- 

Example 4.2. This example shows that transitions at different sites are not independent in general. By 
comparing one CPA time step to one time step with the discretized FPO we see that a specific local 
transition at one site cannot take place if another specific local transition happens at a neighboring site. 
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Figure 4: Illustration of a transition from x = (4, 4, 2, 2) to = (2, 2, 0, 2) in Ex. 14.21 The left side shows 
the preimage, the right the image state. The horizontal numbers correspond to the respective sites, while 
the vertical numbers display E ~ {0, ...,4}. The states x ^-nd V' a-re marked by black rectangles at the 
corresponding sites. 



although both transitions are allowed locally. Consider / = {!,..., 4}, and the system on dynamically 

invariant state space fi™ = [0, 1]"' given for all g N by u"^^ = , = ^^"^^ for i e {1, m— 1}. 
We have U = {0, 1} and define 5 intervals 

rij = [wj,u)j+i), for j G {0, 3}, Vti = [wi,'wz] 

with 

wq = 0, wi ~ 0.183, W2 = 0.31, W3 = 0.4, ^4 = 0.7, wr-, = 1, 

name them by their index and obtain a partition of 57 with E ~ {0, 4}, see Fig. |4l The induced flow is 
denoted by for one time step. We consider the CPA with = J7 for deterministic input g G Dbc(E^) 
given by g{ip) = 5^^^, where x = (4, 4, 2, 2) S E^ . We focus on the image state ij; = (2, 2, 0, 2) g E^ and 
determine 

^ 0.8, ?2 = 0.3625, ri = 0.83875, = 0.32375 
as the solution of the equations 

h{wi,li) ^ w^i, h{lij2)=W2, h{ri,r2) = W2, ft.(r2, W2) = Wi- 
It is possible to show that 

{^X\i + u+v n $"^(f7v|i + v)} ^{v en\w4<V2 <h,l2 <V3 < W3}, 

{^Mx\2+u+v) n $-\r!ai(V|2+v))} '^{ven\ri<V2<l,W2<V3< rs}, 
./o(CTi(x|i+(7+y))(cri(i/'|i+y))) > and fo{a2{x\2+u+v)){cr2{ij\2+v)) > 0. Hence af/3{g){ip) > 0, but 



A(0) 



0, 



mx) 

and so afP{g) ^ Psig)- 

Both examples are scalable in the sense that we can find analogous partitions of [0,c],c e (0, 1), with 
the above properties by dividing all phase space coordinates by c and complete the partition in [c, 1] 
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arbitrarily. So for decreasing size of the coding domains we can still find a partition of [0, 1] with the 
above effects: locality errors are independent from resolution errors. 

Furthermore the examples suggest to choose large V for good approximations. Accordingly it can be 
proven that for maximal V the CPA exactly corresponds to the discretized FPO. 

Proposition 4.4. For V = Knax we find that d//3 = Pb- 



4.3 Consistency Theorem 

The preceding results can be composed to a consistency result for uncertainty propagation with CPA: 
Up to technical postprocessing time evolution of probability densities with CPA can approximate time 
evolution with the FPO arbitrarily close. 

Corollary 4.5. Let T : fl E he a, uniform partition with resolution Afl and g <E k^^{Dbc{E^))- 
Then 

Proof Because there is only a finite number of possible the limit V — ?> Vmax is taken in the discrete 
topology. So the limit point is the evaluation with Knax, and Prop. 14.41 can be used. In a second step 
Thm. 14.21 can be applied to P^^^g): 

^ aIi^o " v^v^ f^B^^o-ffiKBig) - P<i>-{g)\\i 
= yym^\\RP^.{g)-P^r(g)\\, 
=0. □ 



5 Example 

In this section we comment on how to implement an algorithm to evolve uncertainties with CPA. The 
algorithm is tested at the problem of arsenate transportation and asorption in drinking water pipes, 
and the results are compared to a Monte Carlo calculation. Although the theory has been developed 
for uncertainties in initial conditions, in this applicational part we extend the concept slightly such that 
CPA can cope with certain stochastic boundary conditions. This is important in contaminant transport 
modeling [43]. With this first generalization we want to provide evidence that with CPA the treatment 
of more general stochastic spatio-temporal processes seems feasible. 

5.1 Stochastic Boundary Conditions 

We deal with stationary temporal white noise boundary conditions 

gi e D{E'^') for Ki = {!,..., r}, 

gr e D{E^-) for Kr = {m - s + 1, ra] 

instead of deterministic p G E^ . With stationary we mean that the densities do not change in time, and 
the term temporal white noise indicates that there are no correlations in the boundary random variable's 
realizations at different times. 

For this purpose the global function in Def. 13.91 is extendend to 

f ■■ DdBe ^ DdBe, 9 ^ f (g) , 
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5(0 M = 9i{x): 

x^e'^I B.t. (x)l{i,...,^+i,_i} = 
XEE-Kr a.t. »i _ (x ) I { „ _ 3+1 _ i + „) = 

fi9)ii)W = 51 ■"£7(i)(f^'(5)l(7(»))('/'ly+;7(o) ■5(«)(<^) ' foi^)W- 

Furthermore the relations between global and de Bruijn densities then have to be generalized to a : 
D{E^') X D^as^ X D{E^'-) ^ D{E'), {gi xgx .g,) ^ a{{gi xgx g,)) with 

Him X .9 X 9r)){ll^) = .9i(V'k,)a/(ff)('0l{l+r,...,m-s})3r(V'kJ 

and /3 : D{E') ^ D^E''') x olse x 3 ^ Hg) = (.9;, 3/, 5.) with 

xeE' s. t. xUv=V' 

xeE-f B. t. 
x|i+v-=<T-i(iA) 

9rW = Y six)- 

XdE' s. t. x\ki=-4> 

CPA with stochastic boundary conditions may be used to approximate spatio-temporal processes with 
deterministic dynamics, in which the initial and boundary conditions are stochastic. We remark that it 
is straight-forward to use time-dependent stochastic boundary conditions instead of stationary ones. 

5.2 Implementation 

From an implcmcntational point of view two steps of uncertainty propagation with CPA have to be distin- 
guished. Step one is the tranlation of the completely continuous system into a CPA. This is independent 
of initial or boundary conditions and can be achieved in a preprocessing procedure. Step two consists of 
the CPA evolution with given initial and boundary values. It turns out that step one is numerically more 
expensive than step two. For industrial applications like simulation based system monitoring the CPA 
method points towards real-time uncertainty quantification, because the slow step one only has to be 
performed once before the actual simulation. We furthermore note that by construction the simulation 
is parallclizablc in space. 

Step one basically consists of the approximation of local transition probabilities. We propose a local 
version of the standard Monte Carlo quadrature approach |18| in set oriented numerics for this purpose. 
We remark that also advanced adaptive methods have been suggested, see e. g. |13j . 

i) For if e £'^+^ choose test vectors Wi = {wi-r-p, ■■■,'Wi^s+q) & (R")^+^, where {wij}i<w^ is 
randomly distributed over coding domain ^^p^j) ^ ^, respectively. 

ii) Compute for all i < W^p the image points 

Wi = {h{T, Wi-r-p, Wi^s-p), hlj, Wi^^r-p+l, Wi^s-p+l), h{T, Wi^^r+q, Wi^s+q))- 

iii) Determine ipi, ...jipL G E^ such that there is I < L and Wi with T{{wi)j) = for all j £ V. Let 
the number of image points in the specific coding domain be denoted by , i.e. X^^^i ^4>i — 



18 



The local transition function is then approximated by 



W.^IW^ forall VgWi,-^^} 
else 



With increasing number of test points the approximation is expected to get better. However, the number 
of evaluations grows exponentially in \y\. So we suggest to use dc Bruijn calculus to determine transition 
probabilities for large V from transition probabilities for smaller see Fig. [3l For given /o : E^^^' 
D{E^) and W given hyV^^V + W 

where 

/o(<p)(V') = otw{g){ip), 

geiDiE^yr, m = M^Uv+u)- 
It can be shown with an example similar to Ex. 14.21 that this is again just an approximation of the 
directly calculated /q. 

Regarding step two, the simulation with CPA, we remark that it is important to only follow and store 
states with probability larger than a specified threshold whenever possible. Otherwise already for rea- 
sonably large E oi V the calculations are not feasible. An example is the de Bruijn density D^^^, where 
|£^|l^l numbers would have to be handled at every site in I. The set of states with positive probability is 
much smaller, although it typically first grows and then shrinks again in the transient phase of dynamics. 
Note that our de Bruijn choice of aw enables such sparse calculations, whereas the whole space is needed 
to solve for example a linear nonnegative least squares problem. 

5.3 Arsenate Fate in Water Pipe 

Consider the advection and adsorption of arsenate in drinking water pipes, a topic that has attracted a 
lot of attention in the water supply community lately |38l [25] . We describe a water tank on a hill and a 
pipe to a consumer in a valley. Report locations to observe the arsenate concentrations are installed in 
a distance of Ax, see Fig. [5al The physics is described by the Langmuir adsorption model [27] 

dtD + vd^D = -— 1 ^ 1 ,1 Tz{D{S„^ax -A)- KeqA), 

= 1 ^ ^AD{Srna. - A) ~ K,,A) 

ki + kj K'^m.ax A) 

where D is the concentration of dissolved arsenate in ^ and A the concentration of arsenate adsorbed 
at the pipe wall in ^ . We adopt realistic parameter values from [55] [33] 

m I 

V = 10 — —, Th = 50^, 
mm 

fcl = 0.2 —, S,nax = 100^, 

mg mm 

= 0.0537^, kf = 2.4 ,V , 

t 771^ mm 



and consider the system on the approximately positively invariant f2 given by D S [0, 1] and A € [0, Smax]- 
The backward difference with U = {—1, 0}, Ax = 100m and At = Ax/v = 10mm is used. 
To obtain the local function of a CPA we map test points by using intermediate steps on the basis of 
the method of characteristics with the smaller At' = 0.1mm and Ax' = Im. We use V = V ^ {0} and 
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Figure 5: (a) A reservoir on a hill is connected to a consumer in a valley through a pipe with 6 report 
locations, (b) The phase space at every report location is divided into 5x5 coding domains, and the 
steady states are drawn in green. For example /o(((l,4), (2,4))) can be approximated by the transition 
of blue test points in domain (1,4) and black ones in domain (2,4) to the set of red points, (c) shows 
the initial conditions for an exemplary simulation with the according CPA, and the results after 1 and 7 
hours are shown in (d) and (e), respectively. See also Fig. El 
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Figure 6: Steady states after 24 hours, (a) shows the result of a Monte Carlo computation. In (b) one 
finds the result for a CPA with a state space resolution of 5 x 5 and pattern size V = V = {0}, in (c) the 
results for 5x15 and V = V = {0}, and in (d) the ones for 5 x 5 and V = {0}, V = {-1, 0}. 



partition the phase space equidistantly with 5 symbols in each of the n — 2 directions. If we label the 
coding domains from to 4 in each direction, the corresponding CPA results from transition probabilities 
like 

( 0.806 ifV' = (l,4) 
/o(((l,4),(2,4)))(V') = <^ 0.194 ifV = (2,4) , 
[ else 

see Fig. I5bl White noise boundary conditions are applied to describe a random arsenate source in the 
tank, and deterministic initial values represent a pipe which is completely empty in the beginning. The 
observed dynamics is shown in Fig. [5cll5el Dissolved arsenate is transported along the pipe, and over time 
the walls are covered more and more with adsorbed arsenate. After 24 hours a steady state is reached, 
and we compare it to a Monte Carlo calculation, see Fig. I6all6bl The latter has also been obtained on 
the basis of the method of characteristics with At' = O.lmin and Ax' = Im for 20000 evaluations. The 
boundary condition has been drawn from the stationary boundary distribution every lOmm and held 
constant in the meantime. Our example features an interesting probabilitic effect due to the nonlinearity 
of the reaction equations. Although the boundary values arc distributed in the D-domains 2 — 4, the 
consumer mostly observes dissolved arsenate at a concentration of domain 4 in the steady state. 
Furthermore we plot the steady state results from CPA for which the approximation parameters are 
altered. In Fig. |6c] the result is plotted for = y = {0} with an equidistant phase space partition 
of 5 domains in D- and 15 domains in A-direction, whereas in Fig. I6dl the pattern size is extended by 
W = {—1,0} to V = {—1,0}. It is observed that in this example increasing the pattern size does not 
improve the result if compared to the Monte Carlo case, but increasing the state space resolution has a 
notable effect. In all cases we used 75 test points for the coding domain at site and 37 at site —1 to 
approximate the local function, and a probability threshold of 0.00005 in the simulation. 
Wc note that there is often no interest in global results and accordingly no need to transform between local 
and fully global states with a. Local information like indicated in the graphs can be directly extracted 
from the CPA result. Similarly, in practice the information about initial values is often given locally, 
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such that there is no need to use the fuh (3. Besides, we note that the discrete state space information 
is often completely sufficient in practice. In the example a consumer is interested rather in risk level or 
threshold information about water contamination than in information in form of exact concentrations. 
In some biological sytems even the Boolean case, \E\ ~ 2, is enough [3]. 

6 Conclusion 

We have introduced a numerical scheme for density based uncertainty propagation in certain spatio- 
temporal systems. It consists of a preprocessing step, in which the underlying PDE system is translated 
into a cellular probabilistic automaton, and a simulation step, in which initial and boundary conditions are 
evolved. The simulation is parallelizable in the space extension and fast in relation to the preprocessing. 
Furthermore it computes on discrete states instead of on the continuous phase space. Because the discrete 
states can be interpreted as risk levels, fast uncertainty propagation directly on this simplified state space 
suites industrial demands. 

The algorithm is based on state space discretization like in set oriented numerics and on the de Bruijn state 
idea from cellular automata theory. There are two parameters that allow to control the approximation 
of the exact density evolution: state space resolution and de Bruijn pattern size. We have proven 
consistency of the method for uncertain initial conditions under deterministic dynamics and have paved 
the way towards the handling of spatio-temporal processes with more involved stochastic influence. More 
precisely, it has been shown how to deal with white noise boundary conditions, an important topic for 
example in contaminant transport modeling. However, it seems difficult to preserve temporal correlations 
in random parameters with our algorithm. Future research will focus on white noise parameters. We 
are confident that they only extend the preprocessing, whereas the simulation step is not changed. In 
this sense CPA promise to overcome the curse of dimension in parameter space. Besides, we want to 
investigate how our algorithm performs in more complex applications like the simulation of drinking or 
waste water grids. 
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